knitr::opts_chunk$set(echo = TRUE)
rm(list = ls())
#install.packages("ona", repos = c("https://cran.qe-libs.org", "https://cran.rstudio.org")) # only run once
#install.packages("tma", repos = c("https://cran.qe-libs.org", "https://cran.rstudio.org")) # only run once
most recent data Ayano shared
data <- read.csv("~/Rprojects/ayano/PES_010v6(selected).csv")
units <- c("Speaker", "Role", "TimePhase")
codes_1 <- c("Functions", "Aesthetics", "User", "Vision", "Prototyping") # design actions codes
codes_2 <- c("CA", "ALoK", "CSU", "GCA", "Projective","Regulative","Relational") # shared epistemic agency codes
hoo_rules <- conversation_rules(
(TimePhase %in% UNIT$TimePhase & Role %in% UNIT$Role)
)
notes: accum_1 and set_1 are for design action codes; accum_2 and set_2 are for shared epistemic agency codes
accum_1 <-
contexts(data,
units_by = units,
hoo_rules = hoo_rules) %>%
accumulate_contexts(codes = codes_1,
decay.function = decay(simple_window, window_size = 4),
return.ena.set = FALSE, norm.by = NULL)
Note: Not doing means rotation so I can compare everything in the same space
set_1 <-
# model(accum_1,
# rotate.using = "mean",
# rotation.params =
# list(engineers=accum_1$meta.data$Role=="Engineer",
# servicedesigners=accum_1$meta.data$Role=="ServiceDesigner"))
model(accum_1)
set_1 <-
# model(accum_1,
# rotate.using = "mean",
# rotation.params =
# list(engineers=accum_1$meta.data$Role=="Engineer",
# servicedesigners=accum_1$meta.data$Role=="ServiceDesigner"))
model(accum_1)
set_2 <-
# model(accum_2,
# rotate.using = "mean",
# rotation.params =
# list(engineers=accum_2$meta.data$Role=="Engineer",
# servicedesigners=accum_2$meta.data$Role=="ServiceDesigner"))
model(accum_2)
global visual parameter (let’s make sure all the plots use the same level of multiplier)
set_2 <-
# model(accum_2,
# rotate.using = "mean",
# rotation.params =
# list(engineers=accum_2$meta.data$Role=="Engineer",
# servicedesigners=accum_2$meta.data$Role=="ServiceDesigner"))
model(accum_2)
Model 1: design action codes
Plot 1a: Overall Group Comparisons
traces = c(2:10)
plot(set_1, title = "Groups") |>
units(
points= set_1$points$Role$ProductDesigner,
point_position_multiplier = point_position_multiplier,
points_color = c("red"),
show_mean = TRUE, show_points = FALSE, with_ci = TRUE) |>
units(
points= set_1$points$Role$Engineer,
point_position_multiplier = point_position_multiplier,
points_color = c("blue"),
show_mean = TRUE, show_points = FALSE, with_ci = TRUE) |>
units(
points= set_1$points$Role$ServiceDesigner,
point_position_multiplier = point_position_multiplier,
points_color = c("green"),
show_mean = TRUE, show_points = FALSE, with_ci = TRUE) |>
nodes(node_size_multiplier = 0.3,
node_position_multiplier = node_position_multiplier,
self_connection_color = c("black"))|>
plotly::layout(showlegend = TRUE, legend = list(x = 100, y = 0.9)) |>
plotly::style(name = "Product Designer", traces = traces[1]) |>
plotly::style(name = "Engineer", traces = traces[2]) |>
plotly::style(name = "Service Designer", traces = traces[3])
NA
Plot 1b: Group subtractions
prod_pts = set_1$points$Role$ProductDesigner
eng_pts = set_1$points$Role$Engineer
serv_pts = set_1$points$Role$ServiceDesigner
prod_mean_net = set_1$line.weights %>% dplyr::filter(Role == "ProductDesigner") %>% colMeans()
eng_mean_net = set_1$line.weights %>% dplyr::filter(Role == "Engineer") %>% colMeans()
serv_mean_net = set_1$line.weights %>% dplyr::filter(Role == "ServiceDesigner") %>% colMeans()
plot(set_1, title = "Product Designer vs Engineer") |>
units(
points = prod_pts,
points_color = "red",
show_mean = TRUE, show_points = FALSE, with_ci = TRUE) |>
units(
points = eng_pts,
points_color = "blue",
show_mean = TRUE, show_points = FALSE, with_ci = TRUE) |>
edges(
weights = prod_mean_net - eng_mean_net, # optional multiplier to adjust for readability
edge_size_multiplier = edge_size_multiplier,
edge_arrow_saturation_multiplier = edge_arrow_saturation_multiplier,
node_position_multiplier = node_position_multiplier,
edge_color = c("red","blue")) |>
nodes(
node_size_multiplier = node_size_multiplier,
node_position_multiplier = node_position_multiplier,
self_connection_color = c("red","blue"))
plot(set_1, title = "Product Designer vs ServiceDesigner") |>
units(
points = prod_pts,
points_color = "red",
show_mean = TRUE, show_points = FALSE, with_ci = TRUE) |>
units(
points =serv_pts,
points_color = "green",
show_mean = TRUE, show_points = FALSE, with_ci = TRUE) |>
edges(
weights = prod_mean_net - serv_mean_net, # optional multiplier to adjust for readability
edge_size_multiplier = edge_size_multiplier,
edge_arrow_saturation_multiplier = edge_arrow_saturation_multiplier,
node_position_multiplier = node_position_multiplier,
edge_color = c("red","green")) |>
nodes(
node_size_multiplier = node_size_multiplier,
node_position_multiplier = node_position_multiplier,
self_connection_color = c("red","green"))
plot(set_1, title = "Service Designer vs Engineer") |>
units(
points = serv_pts,
points_color = "green",
show_mean = TRUE, show_points = FALSE, with_ci = TRUE) |>
units(
points = eng_pts,
points_color = "blue",
show_mean = TRUE, show_points = FALSE, with_ci = TRUE) |>
edges(
weights = serv_mean_net - eng_mean_net, # optional multiplier to adjust for readability
edge_size_multiplier = edge_size_multiplier,
edge_arrow_saturation_multiplier = edge_arrow_saturation_multiplier,
node_position_multiplier = node_position_multiplier,
edge_color = c("green","blue")) |>
nodes(
node_size_multiplier = node_size_multiplier,
node_position_multiplier = node_position_multiplier,
self_connection_color = c("green","blue"))
Plot 2: All Means by phase
#get list of time phases and groups loop over or use map function
phases = unique(set_1$points$TimePhase)
groups = unique(set_1$points$Role)
col_list_2 = c("red",#phase 1
"blue", #phase 2
"green") #phase 3
traces = c(2:10)
x = plot(set_1, title = "Group Means by Phase")
count = 1
for (i in 1:length(phases)){
for (j in 1:length(groups)){
points = set_1$points %>% filter(TimePhase == phases[i], Role == groups[j])
x = x |>
units(
points = points,
point_position_multiplier = point_position_multiplier,
points_color = col_list_2[i],
show_mean = TRUE, show_points = F, with_ci = FALSE
) |>
nodes(node_size_multiplier = 0.3,
node_position_multiplier = node_position_multiplier,
self_connection_color = c("black"))|>
plotly::layout(showlegend = TRUE, legend = list(x = 100, y = 0.9)) |>
plotly::style(name = paste0(phases[i]," - ",groups[j]), traces = traces[count])
count = count + 1
}
}
x
NA
Statistical Tests:
In these data, the observations are nested within participants, so we need to check if this nesting has a significant effect. We do this by constructing a confidence interval around the intraclass correlation coefficient (https://en.wikipedia.org/wiki/Intraclass_correlation). If the interval contains zero, nesting is not significant.
reg_dat = set_1$points %>% filter(ENA_DIRECTION == "response")
ICC::ICCest(Speaker,SVD1,reg_dat) #not significant
Warning: 'x' has been coerced to a factor
$ICC
[1] 0.05242903
$LowerCI
[1] -0.2098135
$UpperCI
[1] 0.5007296
$N
[1] 10
$k
[1] 3.272727
$varw
[1] 0.3002277
$vara
[1] 0.01661158
ICC::ICCest(Speaker,SVD2,reg_dat) #not significant
Warning: 'x' has been coerced to a factor
$ICC
[1] -0.1818112
$LowerCI
[1] -0.3333495
$UpperCI
[1] 0.1964787
$N
[1] 10
$k
[1] 3.272727
$varw
[1] 0.07137526
$vara
[1] -0.01098045
Nesting is not significant so we proceed with an OLS regression. Below we set up two models for each dimension and test whether including interaction effects in the analysis significantly improves the modes.
mod_x = lm(SVD1 ~ as.factor(TimePhase) + Role, data = reg_dat)
mod_y = lm(SVD2 ~ as.factor(TimePhase) + Role, data = reg_dat)
mod_x_2 = lm(SVD1 ~ as.factor(TimePhase)*Role, data = reg_dat)
anova(mod_x,mod_x_2)
Analysis of Variance Table
Model 1: SVD1 ~ as.factor(TimePhase) + Role
Model 2: SVD1 ~ as.factor(TimePhase) * Role
Res.Df RSS Df Sum of Sq F Pr(>F)
1 28 2.6341
2 24 1.7559 4 0.87815 3.0006 0.03851 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
mod_y_2 = lm(SVD2 ~ as.factor(TimePhase)*Role, data = reg_dat)
anova(mod_y,mod_y_2)
Analysis of Variance Table
Model 1: SVD2 ~ as.factor(TimePhase) + Role
Model 2: SVD2 ~ as.factor(TimePhase) * Role
Res.Df RSS Df Sum of Sq F Pr(>F)
1 28 1.1117
2 24 0.2156 4 0.89614 24.939 3.02e-08 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Interactions are significant, so we should use the interaction models. Let’s view those models
summary(mod_x_2)
Call:
lm(formula = SVD1 ~ as.factor(TimePhase) * Role, data = reg_dat)
Residuals:
Min 1Q Median 3Q Max
-0.85356 -0.01457 0.00790 0.08624 0.32816
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.1652 0.1352 -1.221 0.23385
as.factor(TimePhase)2 0.7013 0.2066 3.395 0.00239 **
as.factor(TimePhase)3 0.6249 0.2066 3.025 0.00585 **
RoleProductDesigner -0.5044 0.1913 -2.637 0.01443 *
RoleServiceDesigner -0.6041 0.1913 -3.158 0.00425 **
as.factor(TimePhase)2:RoleProductDesigner 0.1937 0.2815 0.688 0.49803
as.factor(TimePhase)3:RoleProductDesigner 0.6236 0.2922 2.134 0.04324 *
as.factor(TimePhase)2:RoleServiceDesigner -0.2347 0.2815 -0.834 0.41272
as.factor(TimePhase)3:RoleServiceDesigner 0.6444 0.2815 2.289 0.03117 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.2705 on 24 degrees of freedom
Multiple R-squared: 0.8261, Adjusted R-squared: 0.7681
F-statistic: 14.25 on 8 and 24 DF, p-value: 2.062e-07
summary(mod_y_2)
Call:
lm(formula = SVD2 ~ as.factor(TimePhase) * Role, data = reg_dat)
Residuals:
Min 1Q Median 3Q Max
-0.196826 -0.035322 -0.001122 0.020930 0.273378
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.17794 0.04739 -3.755 0.000977 ***
as.factor(TimePhase)2 -0.08890 0.07239 -1.228 0.231326
as.factor(TimePhase)3 0.62627 0.07239 8.651 7.68e-09 ***
RoleProductDesigner 0.24264 0.06702 3.620 0.001366 **
RoleServiceDesigner 0.29823 0.06702 4.450 0.000168 ***
as.factor(TimePhase)2:RoleProductDesigner 0.06701 0.09865 0.679 0.503477
as.factor(TimePhase)3:RoleProductDesigner -0.35748 0.10237 -3.492 0.001880 **
as.factor(TimePhase)2:RoleServiceDesigner -0.27608 0.09865 -2.799 0.009962 **
as.factor(TimePhase)3:RoleServiceDesigner -0.93799 0.09865 -9.508 1.30e-09 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.09478 on 24 degrees of freedom
Multiple R-squared: 0.89, Adjusted R-squared: 0.8534
F-statistic: 24.28 on 8 and 24 DF, p-value: 1.032e-09
Viewing interactions
emmip(mod_x_2, Role ~ TimePhase)
emmip(mod_y_2, Role ~ TimePhase)
When interactions are introduced, interpreting regression coefficients becomes more difficult, so we calculated the marginal means and various contrasts (i.e., comparisons between groups of interest) to address our research questions.
RQ1 (X dimension - Design Actions): Here we want to compare the “main effect” of Role on the x dimension—i.e., are there differences between the roles on the X dimension when averaging over TimePhases.
emm_x1 = emmeans(mod_x_2, specs = pairwise ~ Role, weights = "proportional")
NOTE: Results may be misleading due to involvement in interactions
print(emm_x1$contrasts)
contrast estimate SE df t.ratio p.value
Engineer - ProductDesigner 0.251 0.119 24 2.116 0.1077
Engineer - ServiceDesigner 0.487 0.116 24 4.192 0.0009
ProductDesigner - ServiceDesigner 0.236 0.113 24 2.085 0.1142
Results are averaged over the levels of: TimePhase
P value adjustment: tukey method for comparing a family of 3 estimates
Only significant difference on the X is between Engineer and Service designer. Let’s calculate the effect size (difference in standard deviations) of that difference. We will use Cohen’s d.
cohensd = function(diff_,g1,g2){
diff_/(sqrt((sd(g1)^2 + sd(g2)^2)/2))
}
contrasts = as_tibble(emm_x1$contrasts)
diff_ = contrasts$estimate[2]
g1 = reg_dat %>% filter(Role == "Engineer") %>% select(SVD1)
g2 = reg_dat %>% filter(Role == "ServiceDesigner") %>% select(SVD1)
cohensd(diff_ = diff_,g1 = g1$SVD1,g2 = g2$SVD1)
[1] 0.9779281
RQ1 (Y dimension - Design Actions): Here we want to compare the “main effect” of Role on the y dimension—i.e., are there differences between the roles on the X dimension when averaging over TimePhases.
emm_y1 = emmeans(mod_y_2, specs = pairwise ~ Role, weights = "proportional")
NOTE: Results may be misleading due to involvement in interactions
print(emm_y1$contrasts)
contrast estimate SE df t.ratio p.value
Engineer - ProductDesigner -0.157 0.0415 24 -3.771 0.0026
Engineer - ServiceDesigner 0.078 0.0407 24 1.917 0.1556
ProductDesigner - ServiceDesigner 0.235 0.0397 24 5.915 <.0001
Results are averaged over the levels of: TimePhase
P value adjustment: tukey method for comparing a family of 3 estimates
Significant differences between engineer/product designer and productdesigner/servicedesigner. Let’s get the effect sizes:
contrasts = as_tibble(emm_y1$contrasts)
diff_evp = contrasts$estimate[1]
g1 = reg_dat %>% filter(Role == "Engineer") %>% select(SVD1)
g2 = reg_dat %>% filter(Role == "ProductDesigner") %>% select(SVD1)
d_evp = cohensd(diff_ = diff_,g1 = g1$SVD1,g2 = g2$SVD1)
diff_pvs = contrasts$estimate[3]
g1 = reg_dat %>% filter(Role == "ProductDesigner") %>% select(SVD1)
g2 = reg_dat %>% filter(Role == "ServiceDesigner") %>% select(SVD1)
d_pvs = cohensd(diff_ = diff_,g1 = g1$SVD1,g2 = g2$SVD1)
d_evp
[1] 0.9137988
d_pvs
[1] 0.8072687
Marginal means - x RQ2
Here we use marginal means to test for differences between timephases within groups and groups within timephase and
library(emmeans)
emm_x1a = emmeans(mod_x_2, specs = pairwise ~ TimePhase|Role, weights = "proportional") #timephases w/n groups
emm_x1a = as_tibble(emm_x1a$contrasts)
emm_x2 = emmeans(mod_x_2, specs = pairwise ~ Role|TimePhase, weights = "proportional") #groups w/n timephases
emm_x2 = as_tibble(emm_x2$contrasts)
#see here: https://cran.r-project.org/web/packages/emmeans/vignettes/interactions.html
emm_x1a = emm_x1a %>% filter(p.value < 0.05)
emm_x2 = emm_x2 %>% filter(p.value < 0.05)
emm_x1a
emm_x2
NA
Several contrasts are significant. Let’s calculate the effect sizes
print(list(t1_evp,t1_evs,t2_evs,t2_pvs))
[[1]]
[1] 1.826995
[[2]]
[1] 2.233997
[[3]]
[1] 4.14715
[[4]]
[1] 1.169235
Some of these are too large to be trusted. This is happening because some individuals are highly similar in the data. We will have to acknowledge this as it is likely due to small sample size
Marginal means -Y RQ2: Let’s repeate what we did for X here
emm_y1 = emmeans(mod_y_2, specs = pairwise ~ TimePhase|Role, weights = "proportional")
emm_y1 = as_tibble(emm_y1$contrasts)
emm_y2 = emmeans(mod_y_2, specs = pairwise ~ Role|TimePhase, weights = "proportional")
emm_y2 = as_tibble(emm_y2$contrasts)
emm_y1 = emm_y1 %>% filter(p.value < 0.05)
emm_y2 = emm_y2 %>% filter(p.value < 0.05)
emm_y1
emm_y2
Several contrasts are significant. Let’s calculate the effect sizes
print(list(t1_ep_y,t1_es_y,t2_ep_y,t2_ps_y,t3_es_y,t3_ps_y))
[[1]]
[1] -2.082427
[[2]]
[1] -2.734181
[[3]]
[1] -13.53816
[[4]]
[1] 2.197361
[[5]]
[1] 8.972244
[[6]]
[1] 8.239644
Model diagnostics
library(performance)
library(see)
check_model(mod_x_2)
Model has interaction terms. VIFs might be inflated.
You may check multicollinearity among predictors of a model without interaction terms.
check_model(mod_y_2)
Model has interaction terms. VIFs might be inflated.
You may check multicollinearity among predictors of a model without interaction terms.
Things look ok for the most part
Plot network subtractions (replace TimePhase and Role for different combinations)
mean1 = set_1$points %>% dplyr::filter(TimePhase == 1, Role == "Engineer") #update here
mean2 = set_1$points %>% dplyr::filter(TimePhase == 1, Role == "ServiceDesigner") #update here
mean1_net = set_1$line.weights %>% dplyr::filter(TimePhase == 1, Role == "Engineer") %>% colMeans()
mean2_net = set_1$line.weights %>% dplyr::filter(TimePhase == 1, Role == "ServiceDesigner") %>% colMeans() #update here
plot(set_1, title = "Network subtraction") |>
units(
points = mean1,
points_color = "red",
show_mean = TRUE, show_points = FALSE, with_ci = TRUE) |>
units(
points = mean2,
points_color = "blue",
show_mean = TRUE, show_points = FALSE, with_ci = TRUE) |>
edges(
weights = mean1_net - mean2_net, # optional multiplier to adjust for readability
edge_size_multiplier = edge_size_multiplier,
edge_arrow_saturation_multiplier = edge_arrow_saturation_multiplier,
node_position_multiplier = node_position_multiplier,
edge_color = c("red","blue")) |>
nodes(
node_size_multiplier = node_size_multiplier,
node_position_multiplier = node_position_multiplier,
self_connection_color = c("red","blue"))
Model 2: Epistemic Agency Codes
Plot 1: Mean Network
node_size_multiplier = 0.3 # scale up or down node sizes
node_position_multiplier = 1 # zoom in or out node positions
point_position_multiplier = 1 # zoom in or out the point positions
edge_arrow_saturation_multiplier = 1.5 # adjust the chevron color lighter or darker
edge_size_multiplier = 2 # scale up or down edge sizes
plot(set_2, title = "Overall Mean--Agency") |>
units(
points= set_2$points,
point_position_multiplier = point_position_multiplier,
points_color = c("black"),
show_mean = TRUE, show_points = TRUE, with_ci = TRUE) |>
edges(
weights =set_2$line.weights,
edge_size_multiplier = edge_size_multiplier,
edge_arrow_saturation_multiplier = edge_arrow_saturation_multiplier,
node_position_multiplier = node_position_multiplier,
edge_color = c("black")) |>
nodes(
node_size_multiplier = node_size_multiplier,
node_position_multiplier = node_position_multiplier,
node_labels = TRUE, # change this to FALSE can remove node labels in case you want to add them back in a nicer font or size for your presentations or publications
self_connection_color = c("black"))
Plot 1b: Overall Group Comparisons
plot(set_2, title = "Groups") |>
units(
points= set_2$points$Role$ProductDesigner,
point_position_multiplier = point_position_multiplier,
points_color = c("red"),
show_mean = TRUE, show_points = FALSE, with_ci = TRUE) |>
units(
points= set_2$points$Role$Engineer,
point_position_multiplier = point_position_multiplier,
points_color = c("blue"),
show_mean = TRUE, show_points = FALSE, with_ci = TRUE) |>
units(
points= set_2$points$Role$ServiceDesigner,
point_position_multiplier = point_position_multiplier,
points_color = c("green"),
show_mean = TRUE, show_points = FALSE, with_ci = TRUE) |>
nodes(node_size_multiplier = 0.3,
node_position_multiplier = node_position_multiplier,
self_connection_color = c("black"))|>
plotly::layout(showlegend = TRUE, legend = list(x = 100, y = 0.9)) |>
plotly::style(name = "Product Designer", traces = traces[1]) |>
plotly::style(name = "Engineer", traces = traces[2]) |>
plotly::style(name = "Service Designer", traces = traces[3])
NA
Plot 1c: Group subtractions
prod_pts = set_2$points$Role$ProductDesigner
eng_pts = set_2$points$Role$Engineer
serv_pts = set_2$points$Role$ServiceDesigner
prod_mean_net = set_2$line.weights %>% dplyr::filter(Role == "ProductDesigner") %>% colMeans()
eng_mean_net = set_2$line.weights %>% dplyr::filter(Role == "Engineer") %>% colMeans()
serv_mean_net = set_2$line.weights %>% dplyr::filter(Role == "ServiceDesigner") %>% colMeans()
plot(set_2, title = "Product Designer vs Engineer") |>
units(
points = prod_pts,
points_color = "red",
show_mean = TRUE, show_points = FALSE, with_ci = TRUE) |>
units(
points = eng_pts,
points_color = "blue",
show_mean = TRUE, show_points = FALSE, with_ci = TRUE) |>
edges(
weights = prod_mean_net - eng_mean_net, # optional multiplier to adjust for readability
edge_size_multiplier = edge_size_multiplier,
edge_arrow_saturation_multiplier = edge_arrow_saturation_multiplier,
node_position_multiplier = node_position_multiplier,
edge_color = c("red","blue")) |>
nodes(
node_size_multiplier = node_size_multiplier,
node_position_multiplier = node_position_multiplier,
self_connection_color = c("red","blue"))
plot(set_2, title = "Product Designer vs ServiceDesigner") |>
units(
points = prod_pts,
points_color = "red",
show_mean = TRUE, show_points = FALSE, with_ci = TRUE) |>
units(
points =serv_pts,
points_color = "green",
show_mean = TRUE, show_points = FALSE, with_ci = TRUE) |>
edges(
weights = prod_mean_net - serv_mean_net, # optional multiplier to adjust for readability
edge_size_multiplier = edge_size_multiplier,
edge_arrow_saturation_multiplier = edge_arrow_saturation_multiplier,
node_position_multiplier = node_position_multiplier,
edge_color = c("red","green")) |>
nodes(
node_size_multiplier = node_size_multiplier,
node_position_multiplier = node_position_multiplier,
self_connection_color = c("red","green"))
plot(set_2, title = "Service Designer vs Engineer") |>
units(
points = serv_pts,
points_color = "green",
show_mean = TRUE, show_points = FALSE, with_ci = TRUE) |>
units(
points = eng_pts,
points_color = "blue",
show_mean = TRUE, show_points = FALSE, with_ci = TRUE) |>
edges(
weights = serv_mean_net - eng_mean_net, # optional multiplier to adjust for readability
edge_size_multiplier = 1,
edge_arrow_saturation_multiplier = edge_arrow_saturation_multiplier,
node_position_multiplier = node_position_multiplier,
edge_color = c("green","blue")) |>
nodes(
node_size_multiplier = node_size_multiplier,
node_position_multiplier = node_position_multiplier,
self_connection_color = c("green","blue"))
Plot 2: All Means
#get list of time phases and groups loop over or use map function
phases = unique(set_2$points$TimePhase)
groups = unique(set_2$points$Role)
col_list_2 = c("red",
"blue",
"green")
traces = c(2:10)
x_agency = plot(set_2, title = "Group Means by Phase -- Agency")
count = 1
for (i in 1:length(phases)){
for (j in 1:length(groups)){
points = set_2$points %>% filter(TimePhase == phases[i], Role == groups[j])
#print(head(points))
x_agency = x_agency |>
units(
points = points,
point_position_multiplier = point_position_multiplier,
points_color = col_list_2[i],
show_mean = TRUE, show_points = F, with_ci = FALSE
) |>
nodes(node_size_multiplier = 0.3,
node_position_multiplier = node_position_multiplier,
self_connection_color = c("black"))|>
plotly::layout(showlegend = TRUE, legend = list(x = 100, y = 0.9)) |>
plotly::style(name = paste0(phases[i]," - ",groups[j]), traces = traces[count])
count = count + 1
}
}
x_agency
Check nesting
reg_dat_agency = set_2$points %>% filter(ENA_DIRECTION == "response")
ICC::ICCest(Speaker,SVD1,reg_dat_agency) #not significant
Warning: 'x' has been coerced to a factor
$ICC
[1] 0.02333382
$LowerCI
[1] -0.2268987
$UpperCI
[1] 0.4706009
$N
[1] 10
$k
[1] 3.272727
$varw
[1] 0.2901337
$vara
[1] 0.00693167
ICC::ICCest(Speaker,SVD2,reg_dat_agency) #significant we will need to include a variable for individual in the model
Warning: 'x' has been coerced to a factor
$ICC
[1] 0.5156316
$LowerCI
[1] 0.1639302
$UpperCI
[1] 0.8233863
$N
[1] 10
$k
[1] 3.272727
$varw
[1] 0.08844428
$vara
[1] 0.09415284
mod_x_agency = lm(SVD1 ~ as.factor(TimePhase) + Role, data = reg_dat_agency)
mod_y_agency = lm(SVD2 ~ as.factor(TimePhase) + Role, data = reg_dat_agency)
summary(mod_x_agency)
Call:
lm(formula = SVD1 ~ as.factor(TimePhase) + Role, data = reg_dat_agency)
Residuals:
Min 1Q Median 3Q Max
-0.72619 -0.10243 0.00017 0.18426 0.34417
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.3092 0.1103 -2.804 0.00906 **
as.factor(TimePhase)2 0.6862 0.1201 5.713 3.97e-06 ***
as.factor(TimePhase)3 1.0597 0.1232 8.601 2.40e-09 ***
RoleProductDesigner -0.2519 0.1257 -2.004 0.05486 .
RoleServiceDesigner -0.4308 0.1232 -3.496 0.00159 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.2873 on 28 degrees of freedom
Multiple R-squared: 0.7564, Adjusted R-squared: 0.7216
F-statistic: 21.74 on 4 and 28 DF, p-value: 2.999e-08
summary(mod_y_agency)
Call:
lm(formula = SVD2 ~ as.factor(TimePhase) + Role, data = reg_dat_agency)
Residuals:
Min 1Q Median 3Q Max
-0.4179 -0.2658 -0.1107 0.2753 0.5530
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.12309 0.12362 0.996 0.32793
as.factor(TimePhase)2 0.01307 0.13464 0.097 0.92338
as.factor(TimePhase)3 -0.05497 0.13812 -0.398 0.69364
RoleProductDesigner -0.50503 0.14093 -3.584 0.00127 **
RoleServiceDesigner 0.15828 0.13812 1.146 0.26151
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.3221 on 28 degrees of freedom
Multiple R-squared: 0.4817, Adjusted R-squared: 0.4076
F-statistic: 6.505 on 4 and 28 DF, p-value: 0.0007821
check for interactions: (interactions are significant so we should use those models)
mod_x_agency_2 = lm(SVD1 ~ as.factor(TimePhase) + Role + as.factor(TimePhase)*Role, data = reg_dat_agency)
mod_y_agency_2 = lm(SVD2 ~ as.factor(TimePhase) + Role + as.factor(TimePhase)*Role, data = reg_dat_agency)
summary(mod_x_agency_2)
Call:
lm(formula = SVD1 ~ as.factor(TimePhase) + Role + as.factor(TimePhase) *
Role, data = reg_dat_agency)
Residuals:
Min 1Q Median 3Q Max
-0.72200 -0.04168 0.00319 0.12198 0.33028
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.46990 0.12236 -3.840 0.000788 ***
as.factor(TimePhase)2 1.18141 0.18691 6.321 1.56e-06 ***
as.factor(TimePhase)3 1.10001 0.18691 5.885 4.52e-06 ***
RoleProductDesigner -0.03821 0.17304 -0.221 0.827107
RoleServiceDesigner -0.16249 0.17304 -0.939 0.357071
as.factor(TimePhase)2:RoleProductDesigner -0.55246 0.25471 -2.169 0.040216 *
as.factor(TimePhase)3:RoleProductDesigner -0.15842 0.26433 -0.599 0.554563
as.factor(TimePhase)2:RoleServiceDesigner -0.84961 0.25471 -3.336 0.002761 **
as.factor(TimePhase)3:RoleServiceDesigner -0.00878 0.25471 -0.034 0.972787
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.2447 on 24 degrees of freedom
Multiple R-squared: 0.8485, Adjusted R-squared: 0.798
F-statistic: 16.8 on 8 and 24 DF, p-value: 4.228e-08
summary(mod_y_agency_2)
Call:
lm(formula = SVD2 ~ as.factor(TimePhase) + Role + as.factor(TimePhase) *
Role, data = reg_dat_agency)
Residuals:
Min 1Q Median 3Q Max
-0.39436 -0.08520 0.00161 0.07954 0.63608
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.09125 0.11044 -0.826 0.416831
as.factor(TimePhase)2 0.58628 0.16871 3.475 0.001959 **
as.factor(TimePhase)3 0.08627 0.16871 0.511 0.613779
RoleProductDesigner 0.03593 0.15619 0.230 0.820034
RoleServiceDesigner 0.26034 0.15619 1.667 0.108555
as.factor(TimePhase)2:RoleProductDesigner -1.01499 0.22991 -4.415 0.000184 ***
as.factor(TimePhase)3:RoleProductDesigner -0.74981 0.23859 -3.143 0.004411 **
as.factor(TimePhase)2:RoleServiceDesigner -0.61494 0.22991 -2.675 0.013251 *
as.factor(TimePhase)3:RoleServiceDesigner 0.23732 0.22991 1.032 0.312246
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.2209 on 24 degrees of freedom
Multiple R-squared: 0.791, Adjusted R-squared: 0.7214
F-statistic: 11.36 on 8 and 24 DF, p-value: 1.662e-06
anova(mod_x_agency,mod_x_agency_2)
Analysis of Variance Table
Model 1: SVD1 ~ as.factor(TimePhase) + Role
Model 2: SVD1 ~ as.factor(TimePhase) + Role + as.factor(TimePhase) * Role
Res.Df RSS Df Sum of Sq F Pr(>F)
1 28 2.3112
2 24 1.4373 4 0.87388 3.648 0.01853 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
anova(mod_y_agency,mod_y_agency_2)
Analysis of Variance Table
Model 1: SVD2 ~ as.factor(TimePhase) + Role
Model 2: SVD2 ~ as.factor(TimePhase) + Role + as.factor(TimePhase) * Role
Res.Df RSS Df Sum of Sq F Pr(>F)
1 28 2.9044
2 24 1.1710 4 1.7333 8.8813 0.0001506 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Need to run mixed model for Y ##model is singular so stick with OLS
library(lmerTest)
#mod_y_agency_mixed = lmerTest::lmer(SVD2 ~ 1 + (1|Speaker) + as.factor(TimePhase)*Role ,data = reg_dat_agency)
#summary(mod_y_agency_mixed)
Marginal means x - RQ1
emm_ag1 = emmeans(mod_x_agency_2, specs = pairwise ~ Role, weights = "proportional")
NOTE: Results may be misleading due to involvement in interactions
print(emm_ag1$emmeans)
Role emmean SE df lower.CL upper.CL
Engineer 0.2572 0.0777 24 0.097 0.4175
ProductDesigner -0.0131 0.0740 24 -0.166 0.1396
ServiceDesigner -0.1911 0.0708 24 -0.337 -0.0449
Results are averaged over the levels of: TimePhase
Confidence level used: 0.95
print(emm_ag1$contrasts)
contrast estimate SE df t.ratio p.value
Engineer - ProductDesigner 0.270 0.107 24 2.520 0.0476
Engineer - ServiceDesigner 0.448 0.105 24 4.265 0.0008
ProductDesigner - ServiceDesigner 0.178 0.102 24 1.737 0.2123
Results are averaged over the levels of: TimePhase
P value adjustment: tukey method for comparing a family of 3 estimates
Viewing interactions
emmip(mod_x_agency_2, Role ~ TimePhase)
Marginal means - x RQ2
emm_ag1 = emmeans(mod_x_agency_2, specs = pairwise ~ TimePhase|Role, weights = "proportional")
emm_ag1 = as_tibble(emm_ag1$contrasts)
emm_ag2 = emmeans(mod_x_agency_2, specs = pairwise ~ Role|TimePhase, weights = "proportional")
emm_ag2 = as_tibble(emm_ag2$contrasts)
emm_ag1 = emm_ag1 %>% filter(p.value < 0.05)
emm_ag2 = emm_ag2 %>% filter(p.value < 0.05)
emm_ag1
emm_ag2
NA
Effect sizes
#todo
Marginal means y - RQ1
emm_ag_y1 = emmeans(mod_y_agency_2, specs = pairwise ~ Role, weights = "proportional")
NOTE: Results may be misleading due to involvement in interactions
print(emm_ag_y1$contrasts)
contrast estimate SE df t.ratio p.value
Engineer - ProductDesigner 0.530 0.0968 24 5.470 <.0001
Engineer - ServiceDesigner -0.127 0.0949 24 -1.341 0.3868
ProductDesigner - ServiceDesigner -0.657 0.0925 24 -7.104 <.0001
Results are averaged over the levels of: TimePhase
P value adjustment: tukey method for comparing a family of 3 estimates
Viewing interactions
emmip(mod_y_agency_2, Role ~ TimePhase)
Marginal means - y RQ2
emm_ag_y1 = emmeans(mod_y_agency_2, specs = pairwise ~ TimePhase|Role, weights = "proportional")
emm_ag_y1 = as_tibble(emm_ag_y1$contrasts)
emm_ag_y2 = emmeans(mod_y_agency_2, specs = pairwise ~ Role|TimePhase, weights = "proportional")
emm_ag_y2 = as_tibble(emm_ag_y2$contrasts)
emm_ag_y1 = emm_ag_y1 %>% filter(p.value < 0.05)
emm_ag_y2 = emm_ag_y2 %>% filter(p.value < 0.05)
emm_ag_y1
emm_ag_y2
Calculate eff sizes
#to do
Model diagnostics
check_model(mod_x_agency_2)
Model has interaction terms. VIFs might be inflated.
You may check multicollinearity among predictors of a model without interaction terms.
check_model(mod_y_agency_2)
Model has interaction terms. VIFs might be inflated.
You may check multicollinearity among predictors of a model without interaction terms.
Looks ok
Plot network subtractions
mean1 = set_2$points %>% dplyr::filter(TimePhase == 1, Role == "Engineer")
mean2 = set_2$points %>% dplyr::filter(TimePhase == 1, Role == "ProductDesigner")
mean1_net = set_2$line.weights %>% dplyr::filter(TimePhase == 1, Role == "Engineer") %>% colMeans()
mean2_net = set_2$line.weights %>% dplyr::filter(TimePhase == 1, Role == "ProductDesigner") %>% colMeans()
plot(set_2, title = "Network subtraction") |>
units(
points = mean1,
points_color = "red",
show_mean = TRUE, show_points = FALSE, with_ci = TRUE) |>
units(
points = mean2,
points_color = "blue",
show_mean = TRUE, show_points = FALSE, with_ci = TRUE) |>
edges(
weights = mean1_net - mean2_net, # optional multiplier to adjust for readability
edge_size_multiplier = edge_size_multiplier,
edge_arrow_saturation_multiplier = edge_arrow_saturation_multiplier,
node_position_multiplier = node_position_multiplier,
edge_color = c("red","blue")) |>
nodes(
node_size_multiplier = node_size_multiplier,
node_position_multiplier = node_position_multiplier,
self_connection_color = c("red","blue"))